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Abstract. A Particle-in-Cell (PIC) numerical simulation of the electron Weibel instability is applied in a 
frame of Darwin (radiationless) approximation of the self-consistent fields of sparse plasma. As a result, we 
were able to supplement the classical picture of the instability and, in particular, to obtain the dependency of 
the basic characteristics (the time of development and the maximum field energy) of the thermal anisotropy 
parameter, to trace the dynamic restructuring of current filaments accompanying the nonlinear stage of the 
instability and to trace in detail the evolution of the initial anisotropy of the electron component of plasma. 

Keywords: Weibel instability, Particle-in-Cell, PIC, Vlasov-Darwin model. 
PACS: 52.35.Qz. 



Introduction 

As it is well known, anisotropic distribution of ther- 
mal velocities in homogeneous collisionless plasma may 
lead to Weibel instability (WI) in case of spontaneous 
generation of transverse to direction of anisotropy mag- 
netic field [Tj. The peculiarity of WI is to have a 
regime where unstable growth of magnetic field gen- 
erated by arisen current layers (filaments in multidi- 
mensional case) is followed by dynamic transformation 
of the current structures and stabilization of the mag- 
netic field energy after saturation of the instability. 

The general picture of Weibel instability mentioned 
above has significant variations in several of practically 
important applications of plasma physics, e.g. WI in 
relativistic plasma which absorbed high-energy EM im- 
pulse [2] is considerably different comparing to WI oc- 
curring in the current sheet of Earth's magnetotail [3] , 
which in turn differs from WI of intense ion beams 
[1]. Recent researches indicate that WI might have an 
important role in dense quantum plasmas [5], where 
WI is supposed to be responsible for the generation of 
non-stationary magnetic fields in compact astrophysi- 
cal objects as well as in laser fusion experiments. 

That variety of forms along with fundamental na- 
ture of WI stimulate its further study using different 
approaches, one of which is a PIC-method based on 



Vlasov-Darwin model [6] applied below. Presence of 
well-developed linear theory of electromagnetic insta- 
bilities [7J, a typical representative of which is WI, al- 
lows one to verify the obtained numerical results. 

This work is devoted to the numerical simulation 
of the electron Weibel instability in order to clarify 
its general picture, in particular, the evolution of the 
initial anisotropy on the time, and the dependency of 
key instability characteristics of the initial anisotropy. 
Ions due to inertia are considered to be a neutralizing 
positive background. 

1 Elements of linear theory 
of Weibel instability 

In order to get notion of the generation mechanism 
of the Weibel instability we will review the basic results 
of linear theory applied to it in the most simple, 1D3V 
(x,v x ,v y ,v x ), setup. 

Let's consider a homogeneous collisionless plasma 
with a positive (ion) background and anisotropic dis- 
tribution of thermal velocities of electrons, given by the 
following distribution function: 
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where no is the density of electrons, u x ,u y ,u z are 
the thermal velocities of the electrons along the corre- 
sponding axes, where we assume (without loss of gen- 
erality) U Z > U X — Uy. 

In this case the plasma, obviously, has an excess of 
fast particles moving along z-axis. However, due to the 
symmetry of the distribution, the total current is zero. 

Now assume that there exists spontaneously gener- 
ated (due to thermal fluctuations) non-zero magnetic 
field B = e v B y Q sin (k x x). Then the Lorentz force 
Fl = —evxBo/c will change the trajectory of a parti- 
cle moving along z, which will lead to formation of spa- 
tially separated current sheets. Localization of these 
sheets will coincide with the zeros of the magnetic field 
which caused the deflection of the particles. Hence, we 
should expect formation of two current sheets (with 
opposite signs) per each wavelength of the magnetic 
field perturbation, amplifying the initial magnetic field. 
The growth of the field, which in turn increases den- 
sity of the current sheets, will continue until the major 
part of the particles become significantly magnetized. 
Fig. [I] shows a qualitative illustration of this process. 




Fig. 1. Qualitative representation of the linear phase 
of Weibel instability. Development of current sheets J 
around the zeros of B y (shown as a surface). B x , B z 
considered zero for simplicity. 

In the described setup the perturbed quantities 
are E z (x,t), B y (x,t) (unperturbed values of which are 
zero) and the electron distribution function /(cc, v, t) = 
fo(v) + fi(x,v,t). 

Let the perturbations (as it's common for lin- 
ear analysis in general) to be represented as 
exp (i (k x x — u)t)). Then, substituting the correspond- 
ing values into the Vlasov equation and the field equa- 
tions, we obtain (as it is shown in [8J) a dispersion 
equation 

where w pe = \J Annoe 2 jm is the electron plasma fre- 



quency, A = (u 2 z /u x — 1) is the anisotropy parameter 
and Z is the plasma dispersion function [S]: 
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where £ in general is a complex value. 

The unstable roots of Eq. ([2| (solutions with purely 
imaginary positive to, the standing waves) are in the 
range from kc/uj pe = to kc/io pe = y/~A for any A > 0. 
Thus in the case when the dimensionless wave number 
k x falls into this interval, there is a sharp increase in 
the amplitude of the perturbation of the corresponding 
wavelength, which is a key feature of the initial stage 
of WI. The growth rate 7 = iio has a single maximum 
in the above-mentioned range, which allows to deter- 
mine the most unstable mode. Note that it's hard (if at 
all possible) to get analytical solution of this equation, 
but it is rather easy to solve it numerically with any 
standard method [TU]. The foregoing is illustrated in 
Fig. [2] where the graphs of the growth rate are shown 
for various values of anisotropy parameter with fixed 
accentuated (u z ) thermal velocity. 
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Fig. 2. Growth rate of WI for various values of 
anisotropy at constant accentuated thermal velocity 
u z = Q.l[c}. 

2 Numerical simulations 

Here we will briefly review the basic features of the 
model we chose, as well as it's geometrical setup and 
the main parameters. 

The general physical properties of the system were 
already described above in the section [l] Now we will 
consider 2.5-dimensional setup (x,y,v x ,v y ,v z ), where 
the accentuated component of the thermal velocity 
is taken perpendicularly to the simulation plane x-y. 
Note that in this case the development of the insta- 
bility will now be supported with formation of current 
filaments (beams), rather than layers, as it was in 1.5- 
dimensions. However, the results of the linear theory 
presented above are still valid in this case, if one as- 
sumes k x = k v = k. 



2 



First of all, we should note that low-frequency na- 
ture of the WI itself allows us to use Darwin (radi- 
ationless) approximation of the electromagnetic fields 
which is significantly more efficient from compu- 
tational point of view comparing to more common full 
Maxwell model. 

Thus, in the PIC-simulations performed here, the 
evolution of plasma system inside of the computa- 
tional domain described by the dynamic equations of 
macroparticles (electrons) 

at m p y c J ,q 
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moving in the self-consistent fields determined by the 
Darwin's equations 
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[E = Ei+E t , V x Ei = 0, VE t = 0, 



where E\ and E t are, correspondingly, the longitudi- 
nal (curl-free) and transverse (divergence-free) compo- 
nents of the electric field. 

As it is easy to see, Darwin's approximation differs 
from the full Maxwell description only by the elim- 
inated transverse displacement current, which physi- 
cally means neglecting radiation and transition to fields 
with instantaneous propagation [12 . At the same time 
the system partially keeps inductive effects (associated 
with Faraday's law) and ensures holding of the conti- 
nuity equation due to the presence of longitudinal part 
of the displacement current. 

Then we specify the uniform spatial distribution of 
electrons and singly charged ions, the latter ones are 
considered as a motionless neutralizing background. In 
the expression ([!]) we choose u z = 0.1 [c], u x = u y = 
0.0316 [c]. For these values of thermal velocities the 
anisotropy parameter equals to A = 9. 

The characteristic linear size of the computational 
domain, which we can choose for the simulation, must 
correspond to the wavelength for which the dispersion 
equation (for given values of anisotropy and accentu- 
ated thermal velocity) shows the instability growth rate 
close to maximum (Fig. [2]). 

In particular, solving numerically the equation ^ 
for our reference case Aq — 9 and u z = 0.1 [c], we ob- 
tain the approximate value of 7 max = 0.037 and the 



corresponding value of fc max = 1.2 [w pe /c], in terms of 
the wavelength it will give us A max — 5.2 [c/u> pe \. 

Thus, choosing a rectangular computational do- 
main L x — L y = 25 [c/u;p e ] will allow us to trace in 
detail the development of initial current system and 
its restructuring in nonlinear stage of WI. 

Determination of the total number of particles N p 
is preconditioned by a natural trade-off between hold- 
ing colisionless property of the model in the time- frame 
sufficient for development of the instability and com- 
putational cost of a singe run. Using theoretical esti- 
mations of the collisionless time for big particles (with 
width W > Xd) [13]) as well as test runs, it was found 
that the number of particles N p of order 10 6 allows 
to have the collisionless period at least 5 times longer 
than it is needed for the development of WI. 

Following the work [5], we choose periodic bound- 
ary conditions for both x and y directions. 

Finally, note that in the series of experiments con- 
sidered below variation of the the initial anisotropy pa- 
rameter Aq has been done by changing of u x and u y 
keeping u z constant in all runs to keep non-relativistic 
regime of the simulations. 

3 Discussion of results 

The instability development process is clearly visi- 
ble on the plots of the average over space magnetic field 
energy density versus time for different initial values of 
the anisotropy shown in Fig. [3] 
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Fig. 3. Average magnetic field energy density versus 
time for different values of initial anisotropy (Aq = 
1..A0) at constant accentuated thermal velocity. Com- 
putational domain L x = L y = 25 [c/ui pe \; u z = 0.1 [c]; 
mesh size 128 x 128; 500 particles of each specie per 
cell; time step r = 0.25 [l/w pe ]. 

The initial value of the average energy density is close 
to zero, but by the time t = 100 [l/w pe ], corresponding 
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to appearing of distinguishable areas of current local- 
ization, it is observed a noticeable increase, which is 
the greater the larger is the initial anisotropy of the 
plasma. As noted above, the linear stage ends when 
the system formed a pronounced current structure con- 
sisting of a system of current beams (Fig.[6j time slice 
t = 150 [l/a;p e ]), and when the particles become signif- 
icantly magnetized in average (mean Larmor radius is 
of the order of the current beam radius) . 

Further development of the WI is nonlinear and is 
accompanied by dynamic merging of equidirected cur- 
rent filaments into larger-scale structures, that can be 
seen from a series of time slices shown in Fig. [6] Along 
with that the magnitude of the magnetic field energy 
stabilizes. Here we would note that this picture is quite 
similar to one considered in [14] . where the authors 
study filament recombination which take place on the 
nonlinear stage of WI arising with propagation of rel- 
ativistic electron beam in background plasma (a the- 
ory of relativistic coalescence has been developed in 

ESE5I). 

The merging continues until the velocity distribu- 
tion of the electrons does not become an equilibrium. 
In that context it is interesting to trace the evolution 
of the initial anisotropy of the electron component dur- 
ing the development of WI, including at the non-linear 
stage. The dependency of the anisotropy on the time 
is shown in Fig. [4] From the figure one can see that the 
process of growth and further stabilization of current 
density and magnetic field energy density is accom- 
panied by isotropization (decrease in the anisotropy) 
of the medium, caused by collective processes of for- 
mation and restructuring of current filaments together 
with competitive development of the unstable modes 
of the magnetic field. 
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Fig. 4- Evolution of the anisotropy for different initial 
An. The same setup as for Fig.^ 



It is interesting to note here an non-obvious fact 
that regardless of the initial value of the the anisotropy, 
its value at the instability saturation stage tends to a 
nonzero threshold value (see Fig.[4j. 

The origin for that residual anisotropy is in the 
problem setup, and it is connected with an actual 
finiteness of the perturbation spectrum of any system 
with periodic boundary conditions, due to the fact that 
in such a system waves longer than linear size of the 
domain can not exist. 

This conclusion is confirmed by the results of [17j . 
where it is shown that for a limited system (in the sense 
stated above) the minimum level of the anisotropy at 
the saturation stage of WI can be defined by the fol- 
lowing expression: 

Axnin (^minC/Wpe) • 

In our case the minimal wave number allowed in the 
system is fc m inc/w pe = 2ir/L x , y = 2ir/25 = 0.25, which 
gives A m ; n = 0.063. Some discrepancy of the experi- 
mentally observed levels of residual anisotropy (Fig. [4]) 
and its theoretically predicted value can be explained 
by the fact, that at time t = 500 [l/w pe ] the dominant 
mode of the instability has not yet reached the limit 
determined by the period of the system. As stated 
above, one of the objectives of these experiments was 
to clarify the dependency of the basic characteristics of 
the instability — its growth rate, the characteristic time 
and the maximum energy density of the magnetic field, 
of the initial anisotropy of plasma. 

In this connection let's return to Fig.[3j from which 
one can see the expected effect: the increase of the 
maximum field energy density and instability growth 
rate with the increase of initial anisotropy. Less trivial 
is the profile of this dependency shown in Fig. [5] fast 
increase of the maximum energy density on the inter- 
val A < 10, followed by an equally fast drop-off with 
A > 10. 

Moreover the value of Ao around 25 can be consid- 
ered as an upper threshold of instability in the sense 
that further increase in the initial anisotropy of elec- 
tron component has virtually no effect on the basic 
characteristics of the WI — the maximum value of mag- 
netic field energy density (in case of fixed value of ac- 
centuated thermal velocity u z ). The observed effect 
and its quantitative characteristics are confirmed by 
analytical expression which can be obtained on the ba- 
sis of the expression derived in the work [17] : 
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(5) 



where T z is the temperature of accentuated electron 
component along z-axis. 
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In the just mentioned work the expression ^ has 
been used as an estimation of evolution of the magnetic 
field energy density, where the values of A were being 
taken from the numerical simulation. 

However, in our context, it can be interpreted as an 
analytical dependency of the magnetic field energy den- 
sity (in units of UqT z ) of the current value of anisotropy 
A with its initial value given as Aq. Considering ^ 
that way, its easy to get an interesting for us estima- 
tion of the maximum energy of WI (or, equivalcntly, 
the maximum proportion of the energy of accentu- 
ated component which can be converted into field en- 
ergy) and also the explanation of the existence of upper 
threshold of WI on Ao mentioned above. Indeed, dif- 
ferentiating the right hand side of ^ with respect to 
A and equating the derivative to zero we obtain: 

A -A{2 + Aq) 
3(A + 1)2(^ + 1) U - 

Where from we get the value of the anisotropy coef- 
ficient which delivers the maximum to the expression 
([5]) with fixed value of A : 

A t Ao 
cxtr ~ 2 + A Q ■ 

Substituting the resulting value of A cxt r back into ([5]), 
we find the desired dependence of the maximum mag- 
netic field energy on the value of A : 



A? 
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n T z0 Al _ 
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where w^q 1 is the ^-component of the average over do- 
main initial kinetic energy density of electrons. De- 
pendency ^ together with the experimental values is 
shown in Fig. [5] 




Fig. 5. The dependency of the maximum average mag- 
netic field energy density of the initial anisotropy at 
constant accentuated thermal velocity u z = 0.1 [c]. The 
same setup as for Fig.uA 



Analysis of the last expression gives a number of in- 
teresting provisions which are in good agreement with 
computer simulation experiments. 

First, WBmax has a limit in case of fixed values of 
^0) T z q (u z ) and Aq — > 00 (the last one actually means 
u x ,u y -t 0) 
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corresponding to the maximum possible energy den- 
sity of WI. Particularly, in our case, where uqT z o — 

0. 005 [n e mc 2 ]), u^max ec l uals to 4- 1 ' 10~ 4 [n e mc 2 ], 
which, as it can be seen in Fig. [5j is in good agree- 
ment with the numerical data. 

Second, observed in experiments, the upper thresh- 
old of WI with respect to A is determined by the value 
of Ao, starting from which the fraction j4q/(1 + Ao) 2 
becomes close to unity. It is easy to see that the value 
of Ao close to 25 meets this requirement (indeed, if 
Ao = 25 the fraction equals to 0.925, and if Ao = 50, 

1. e. increased by factor of two, it only slightly increases 
up to 0.96). 

Third, the portion of kinetic energy stored in the ac- 
centuated component of the system which goes for de- 
velopment of Weibel instability depends on the degree 
of the initial anisotropy (Ao) and it can not be higher 
than 1/6. For instance, in our experiments it varies 
from 12.4% (for A = 9) down to 4.4% (for A = 2). 
These data are consistent both with the analytical pre- 
dictions obtained from the expression (j6| and with the 
upper limit of the converted anisotropy energy (around 
10%) obtained in [jjj, where 2D simulation has been 
performed in the frame of full electromagnetic model 
of plasma. 

If we define the characteristic time of the instabil- 
ity (Tj) as the time needed to reach the maximum of 
magnetic field energy density, then from the same Fig. 
[3] one can easily estimate both the value of Tj as well 
as it's dependence of the initial anisotropy Ao- 

Finally, we note good agreement between the nu- 
merical value of the dimcnsionlcss instability growth 
rate (0.034) with the one predicted by the linear theory 
(0.037) for the setup considered above where Ao = 9 
and u z =0.1 [c). 

Conclusion 

Thus, the performed computer experiments allowed 
to supplement the classical picture of the Weibel insta- 
bility. Particularly, we obtained the dependency of the 
characteristic time of its development and the maxi- 
mum magnetic field energy density of the initial value 
of the anisotropy parameter. Also we have traced the 
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dynamic restructuring of the current filaments accom- 
panying nonlinear stage of instability saturation, as 
well as the evolution of the initial anisotropy of the 
electron component of the plasma. 

As a further work it is planned to perform a numer- 
ical study of the dynamics and generation mechanisms 
of two-species Weibel instability, which implies a ki- 
netic representation not only of the electron but also 
the ion plasma component on the time scale by orders 
of magnitude grater than it is for a single-species WI. 
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J z [n e q e c], B xy [m,,uj pe /q^, t = 0[l/u pe ] J z [n e q e c], B xy [rn e uJ pe /q e ], t = 200[l/wpe] 




Fig. 6. The current density J z and the magnetic field B xy at different time moments for 2. 5- dimensional setup 
of WI in case Aq = 9. Computational domain L x — L y — 25 [c/u) pe \; u z =0.1 [c]; mesh size 256 x 256; 1000 
particles of each specie per cell; time step t = 0.25 [l/w pe ]. 
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